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Abstract. When a cluster or nanodroplet bears charge, its structure and 
thermodynamics are altered and, if the charge exceeds a certain limit, the system 
becomes unstable with respect to fragmentation. Some of the key results in this area 
were derived by Raylcigh in the nineteenth century using a continuum model of liquid 
droplets. Here we revisit the topic using a simple particle-based description, presenting 
a systematic case study of how charge affects the physical properties of a Lennard- Jones 
cluster composed of 309 particles. We find that the ability of the cluster to sustain 
charge depends on the number of particles over which the charge is distributed — a 
parameter not included in Rayleigh's analysis. Furthermore, the cluster may fragment 
before the charge is strong enough to drive all charged particles to the surface. The 
charged particles in stable clusters are therefore likely to reside in the cluster's interior 
even without considering solvation effects. 
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1. Introduction 

Charged clusters and small droplets are important in many areas of physical chemistry, 
ranging from aerosol formation in the atmosphere [1] to electrospray processes in the 
laboratory [2]. The presence of charge on a droplet affects both its interactions with 
other species and the physical properties of the individual droplet. Amongst the most 
important questions relating to physical properties are how much charge a cluster or 
droplet can accommodate without becoming unstable on short timescales, and how the 
break-up of the droplet proceeds as this charge limit is approached and exceeded. 

One area where such considerations are crucial is electrospray ionisation mass 
spectrometry. The electrospray technique produces a stream of charged droplets that 
solvate molecules or complexes to be conveyed into the mass spectrometer for analysis. 
Before the analyte's mass-to-charge ratio is measured, it becomes separated from the 
solvent, and some of the droplet's charge is deposited on it. Various models have 
been proposed to describe the mechanism by which the solvent leaves the analyte 
[3, 4, 5] in different circumstances. Computer simulations with atomistic detail have 
been employed to obtain a more explicit picture of the desolvation process and also to 
study the structure of water droplets containing ions [6, 7, 8] and solvated biomolecules 



On the fundamental question of the stability of pure charged droplets, the 
key theoretical work dates back well over a hundred years to Rayleigh [13], who 
considered arbitrary deformations of a uniformly charged, structureless spherical 
droplet. Deformations away from a spherical shape increase the surface area of the 
droplet, thereby increasing its surface energy, but at the same time decrease the 
electrostatic repulsion by spreading the charge out. At a total charge 



known as the Rayleigh limit, the electrostatic contribution to the energy overwhelms 
the surface energy and the droplet becomes mechanically unstable. In Eq. (1), 7 is 
the surface tension of the liquid, eo is the vacuum permittivity and D is the droplet 
diameter. Rayleigh's formula remains of great practical use to the present day. 

To examine the nature of charge-driven instabilities with particle-level detail, 
simulations of small charged clusters have been carried out by adding Coulomb 
interactions to simple generic potentials like Lennard- Jones [14] and Morse [15]. 
However, these studies have mostly concentrated on total charges where the cluster 
is decisively in the unstable regime and rapidly undergoes fragmentation. 

In this article, we take a closer look at the properties of charged clusters with 
an explicit but generic representation of individual particles, concentrating on the 
structure of clusters as the charge-driven instability is approached from below. We 
treat the clusters at equilibrium, setting aside the dynamical aspects of fragmentation 
for now. The results presented here constitute a case study on a cluster of 309 
Lennard- Jones particles to which Coulomb repulsion has been added. 309 is a "magic 
number" for the Lennard- Jones potential; the lowest-energy structure consists of four 



[9, 10, 11, 12]. 
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complete icosahedral shells around a central atom [16]. The present study represents the 
preliminary steps in a more comprehensive ongoing investigation of charged droplets. 

2. Model and methods 

Our cluster of N = 309 particles is bound by the pairwise potential 



V = 4uJ2 

i<j 




N / \ 12 / \ 6 N 



where u and a are the Lennard- Jones pair well-depth and diameter, respectively, 
is the separation of particles % and j, and is the charge on particle %. By defining a 
reduced charge q* — qij (47re cru) 1 / 2 , along with the usual Lennard- Jones reduced energy 
and lengths V* = V/u and r?- = r^/cr, Eq. (2) may be rewritten in the dimensionless 
form 

JV 
i<j 

We will consider clusters where n of the particles carry a charge q* = q* and the 
remaining N — n particles are uncharged, making the overall charge on the cluster 
Q* = nq* . Although q itself should be an integer multiple of the electronic charge 
e, the reduced charge q* also depends on the Lennard- Jones parameters. We will 
therefore treat q* as a continuously variable effective charge. The fact that the particles 
constituting the cluster are represented explicitly and may each carry a charge of only 
q* or means that there is no need to consider an underlying continuum medium with 
its own conductivity. The charge distribution is completely determined by the function 
V* and a small number of external parameters, such as the temperature. Therefore, 
although our model is coarse-grained, its particulate nature makes it more physically 
detailed than a continuum model. 

We simulate the cluster at thermodynamic equilibrium using canonical Monte Carlo 
(MC) simulations. In addition to standard displacements of individual particles, we 
employ trial moves that attempt to swap the charges on randomly chosen charged- 
uncharged pairs of particles while keeping all positions fixed. These steps are accepted 
according to the usual Metropolis criterion [17]. The swap moves are not intended to 
mimic the dynamical process of charge transfer or conductivity that may occur between 
molecules in certain systems. In particular, the swaps are not restricted to pairs of 
adjacent particles. The purpose of the moves is to enhance the efficiency with which 
the equilibrium distribution of charges is reached and sampled, and the moves achieve 
this by, in effect, decoupling the mobility of the charge from the movement of individual 
particles. The same distribution would eventually be achieved by any set of MC moves 
that obey (detailed) balance. 

We will be concerned with the behaviour of the cluster as a function of the 
magnitude of the single-particle charge q* for a given n and, to some extent, on the 
properties with a given charge as a function of temperature T. To accelerate sampling 
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Figure 1. Heat capacity of the neutral LJ309 cluster, obtained by the multiple 
histogram method. 



in these cases, we have employed replica-exchange MC [18] , using either the temperature 
[19] or the charge as the parameter that varies between replicas. In the latter case, the 
acceptance criterion for exchanging configurations Xi and X2 that are initially found 
in replicas where the single-particle charges are q* and q£, respectively, is 

p acc = min {l, exp [(W(X X ) - W(K 2 )) (qf - qf)/kT] } . 

Here, k is Boltzmann's constant and W(X) = Y^^ 7 "*^ ■, where the sum runs over pairs 
of charged particles. 

Any cluster is thermodynamically unstable with respect to vaporisation unless it is 
confined. We therefore constrain the cluster to lie within a spherical container of radius 
R = 6.5a whose centre follows the centre-of-mass of the cluster [20]. This radius is 
large enough not to inhibit large fluctuations in the shape of the cluster and also allows 
fragmentation to occur. We shall return to the consequences of this particular choice of 
radius when discussing the results. 



3. Stability 

We will consider the effect of adding charge to the Lennard- Jones cluster at reduced 
temperatures T* = kT /u = 0.2 and 0.43. Figure 1 shows the canonical constant- volume 
heat capacity C v of the cluster, derived from temperature replica-exchange simulations 
coupled with multiple histogram reweighting [21, 22]. The taller of the two peaks is due 
to the main melting transition [23], meaning that our two temperatures correspond to 
well within the solid-like state and just within the liquid-like regime, respectively. 

At T* = 0.2, the neutral cluster fluctuates uneventfully about its icosahedral global 
minimum [16]. This qualitative behaviour persists if some charge is added, but for 
sufficiently large charge, a small number of charged particles are expelled, leaving a 
well-defined sub-cluster at the centre of the spherical container. The expelled particles 
typically lie close to the container edge, repelled there by the long range Coulomb 
potential. We may therefore count an individual configuration as dissociated if at least 
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Figure 2. (Colour online) Fraction of configurations in which dissociation is detected 
for LJ309 with n = 100 charged particles at temperatures T* = 0.2 (black circles) and 
0.43 (red squares) as a function of the charge per particle. 

one charged particle lies within a small distance, chosen to be a/2, from the container 
wall. The line marked with circles in Fig. 2 shows that the equilibrium fraction of 
dissociated configurations at equilibrium rises sharply at a particular value of the single- 
particle charge, indicating that electrostatic repulsion has become strong enough to 
overcome the cohesive part of the Lennard- Jones potential. We will designate the total 
charge at which the fraction of dissociated configurations passes through 0.5 by <3ma X ; 
beyond this threshold, more than half the equilibrium configurations are dissociated. 

At the liquid-like temperature T* = 0.43, a similar pattern is observed, as shown by 
the squares in Fig. 2. However, since the cluster has now escaped the basin of attraction 
of its icosahedral global minimum and explores a range of disordered structures, the 
transition from an intact cluster to full dissociation is a broader function of the charge. 
We will nevertheless use the same criterion of Q max to quantify the maximum charge 
that the cluster can sustain. 

Importantly, Q max is a function of n, i.e. the maximum total charge that the cluster 
can sustain depends on the number of particles over which the charge is spread. Figure 3 
shows that Q max increases quite dramatically with n, especially at the lower temperature. 
It is not straightforward to compare the Qr of Eq. (1) with Q max , not least because we 
must work at temperatures well below the triple point of T* = 0.694 [24], where the 
bulk surface tension is not known. However, we note that the influence of n discrete 
charges does not enter into Rayleigh's analysis because in that analysis the droplet is 
treated as being homogeneously charged. The dependence of the maximum charge on n 
is potentially important in the context of small electrospray droplets, where the charge 
of the droplet may be due to a relatively small number of ions. Our thermodynamically 
based Q max suggests that such droplets have a greater tendency to dissociate by the 
emission of charged particles than a hypothetical cluster where the same amount of 
charge is spread out over a larger number of particles. 

Increasing the container radius R decreases the value of Q max that is obtained 
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Figure 3. (Colour online) Maximum total charge Q^&x that the cluster can sustain 
at T* = 0.2 (black circles) and T* = 0.43 (red squares) as a function of the number of 
particles n over which the charge is spread. 



because dissociated charged particles can then go further from the remaining sub-cluster, 
thereby lowering the potential energy of the configuration and raising its statistical 
weight relative to the undissociated cluster. However, the decisive trend of increasing 
Qmax with n persists for all R large enough to permit dissociation. 

If the charge is increased beyond Q max , a further structural change may result. 
It is possible for the cluster to split into two or more large sections, which then repel 
each other as far as the container will allow. This subdivision does not happen if n is 
sufficiently small since the electrostatic energy can then be lowered by the emission of 
individual particles without sacrificing the attractive Lennard- Jones interactions in the 
main part of the cluster. However, as n increases, subdivision occurs at a threshold 
increasingly close to Qmax and, for large container radius, may even preempt significant 
loss of individual charged particles. This observation again points to the importance of 
taking into account not only the total charge borne by a cluster but also the number of 
particles over which the charge is divided. One may anticipate that a dynamic study 
would find increased competition between different instability mechanisms at higher n. 

4. Structure 

An important question concerning charged droplets is where the charged particles reside. 
In the case of ions dissolved in water, highly specific effects come into play in determining 
the affinity of ions for the surface of the bulk liquid itself, according to how a given ion 
influences the interactions between the water molecules [25, 26]. In a droplet and in the 
absence of strong surface affinity effects, one might expect ions to be mutually repelled 
to the droplet surface by electrostatic considerations. However, this is not always found 
to be the case [27, 7]. 

Our model is more straightforward and generic than water, allowing us to consider 
the structure of charged droplets without the rather special effects of hydrogen bonding. 
Energetically, the driving forces in our model are towards the optimisation of as many 
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Lennard- Jones pairwise neighbours as possible and toward the maximisation of distance 
between like charges. The global potential energy minimum structure of the cluster, 
obtained using the basin hopping algorithm [28] is always a slightly relaxed version of 
the neutral icosahedron [16]. As n is increased from small values, the charged particles 
preferentially occupy vertices, then edges and finally faces of the outer shell before 
populating the cluster's interior. Counteracting these energetic tendencies at finite 
temperature is the entropy of permuting charged and neutral particles. 

To determine how the charges are divided between the surface and the interior of the 
cluster, we need a method for identifying surface particles. At solid-like temperatures, 
where the cluster retains a clear structure of concentric shells, the outermost shell of 
162 particles can unambiguously be designated as the surface. For liquid-like structures, 
the surface is less clear-cut and its curvature adds to the difficulty of devising a simple 
algorithm to detect the outermost layer. Existing methods to tackle this problem include 
the cone algorithm of Dellago and coworkers [29], which identifies surface particles as 
those for which the vertex of a cone with a particular aperture can be placed at the 
particle's centre without any other particles falling inside the cone. The cone may 
have any orientation, thereby taking into account the fact that the surface normal of 
a droplet may locally deviate significantly from the vector joining that region to the 
droplet 's centre. 

We have devised a somewhat more straightforward algorithm that shares the 
advantageous features of the cone algorithm and produces a reliable and intuitively 
realistic division between surface and core particles. The method borrows an idea from 
algorithms used to determine the surface of proteins [30]. The principle is that a particle 
should be regarded as lying on the surface if it would be accessible by a fictitious probe 
particle approaching from a distance, while the cluster configuration is held fixed. To 
determine this accessibility, a nominal hard-sphere diameter D^j is assigned to the 
Lennard- Jones particles. We then attempt to place a probe particle of diameter D pr 
in contact with any combination of three particles that are sufficiently close for three 
contacts to be achieved simultaneously. For each such triplet, there are two possible 
positions for the probe particle, located symmetrically on either side of the plane defined 
by the triplet on the line passing through the circumcentre of the triangle. If the probe 
particle in at least one of these positions does not overlap with any other particle in the 
cluster then all three Lennard- Jones particles are designated as lying on the surface. A 
given particle may qualify as being on the surface by being tested in conjunction with 
any two of its neighbours. Although a systematic test over all triplets of particles may 
seem computationally expensive, it is quite easy to omit needless iterations of the outer 
loops, and neighbour lists [31] could be implemented for efficiency in larger systems if 
necessary. 

The only parameters to fix in this probe method are the nominal diameters -Dlj and 
_D pr . Both should be comparable with a and the results should not be sensitive to the 
values chosen. We have found that choosing a slightly smaller probe diameter enables 
the algorithm to cope better with crevices in the surface and adopted D LJ = 1.2a, 
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Figure 4. (Colour online) Top: Illustration of the method for distinguishing surface 
particles (light/yellow) from the core of the cluster (dark/blue). The surface particles 
touching the probe (small, red) are highlighted. Bottom: liquid-like snapshot of the 
cluster, highlighting the surface particles. 

.Dpr = O.80" in this study. The probe method and an example of its outcome are 
illustrated in Fig. 4. 

Figure 5 shows how the average number of charged particles occupying surface sites 
varies with the magnitude of the charge for several values of n at T* = 0.43. Due to 
the liquid-like structure, the number of surface sites fluctuates and the sites have to 
be identified for each sample separately. The charge on the horizontal axis of the plot 
has been scaled by the maximum charge that the cluster can sustain for each value 
of n respectively (Fig. 3). Looking at the vertical line Q*/Q^ aa _ x shows that it is only 
for very small numbers of charges and only when the charge approaches the limit of 
stability of the droplet that one finds most of the charged particles on the surface of the 
cluster. Hence, for a cluster of this size bound by the generic Lennard- Jones potential, 
the charge required to drive the majority of particles electrostatically to the cluster 
surface is already too much for the cluster itself to remain intact. Therefore, even in a 
model without specific surface affinity effects or polarisability, charged particles may be 
found predominantly in the interior of the cluster. 
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Figure 5. (Colour online) The fraction of charges that occupy surface sites as a 
function of the total charge of the cluster for different numbers n of charged particles 
at temperature T* = 0.43. For each value of n, the charge axis has been scaled by the 
maximum charge from Fig. 3. The vertical line highlights the limit Q* = Qmax- 



For Q* > Qmax' the loss of charged particles from the cluster causes the number of 
charges at the surface to start decreasing for small n in Fig. 5. However, for higher n, 
although some charged particles are expelled, the higher charge still leads to an overall 
increase in the proportion of charges driven to the surface of the remaining sub-cluster, 
and the upward trend continues beyond Q max - 

In our simple model, the occupation of surface sites by charged particles is 
determined by a competition between electrostatic repulsion of charges to the surface 
and the entropy of mixing of charged and uncharged particles. For solid-like 
temperatures, the balance between these effects can be calculated quite accurately with 
little effort by using a few approximations. At these low temperatures, let us neglect the 
effect of thermal fluctuations on the positions of the particles and consider the structure 
to be the perfect icosahedron of the neutral cluster's global potential energy minimum, 
which has S = 162 surface particles. If n of the N = 309 particles are charged, the 
number of ways of permuting s of them amongst the surface sites and the remaining 
n — s amongst the N — S core sites is 

where C£ = a\/\b\(a — b)\] is a binomial coefficient. Let us also define 

WsM = ('E^- 1 ) > (3) 
\i<j I s 

where the sum runs over the charged particles only and the angle brackets imply an 
unweighted average over configurations in which s charges are randomly assigned to 
surface sites in the frozen cluster and the remaining charges are randomly distributed 
amongst the core sites. Wjv, n ( s ) can be sampled very quickly since the particle positions 
are held fixed while the charges are permuted amongst the sites, and the reciprocal 
pairwise distances need be calculated only once. The mean energy of configurations 
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with s charges at surface sites is then given instantly for any magnitude of the charge 



We may therefore predict the probability that precisely s particles lie on surface sites 
to be 



with the constant of proportionality defined by the normalisation Y%=oPN,n( s ) = 1. 

Figure 6 compares the predictions of Eq. (4) with the results of explicit MC 
simulations. The averages for Wn,u{s) in Eq. (3) were made using 10 5 samples at 
each value of s, which allows Wn jU (s) to be evaluated in a matter of seconds for the 
full range of s. At low n and modest Q*/Q^ 13iX , the predictions of Eq. (4) compare 
essentially perfectly with the simulations even at T* = 0.3, which is approaching the 
first peak in the heat capacity (Fig. 1). As n increases, Eq. (4) begins to overestimate 
the number of charged particles at the surface, but the peak position is only in error 
by a few per cent, even for n = 150 and a total charge approaching Q max . Given that 
Vn n is a crude average and that thermal motion and charge-induced distortion of the 
cluster are neglected in Eq. (4), the simple and fast prediction works well. Importantly, 
we may conclude that the opposing driving forces of electrostatic repulsion and the 
permutational entropy capture most of the effects that determine the location of the 
charges in this model at low temperatures. 

In principle, a similar approach could be taken to the liquid-like state. However, 
the number of surface sites is then no longer fixed and it may be necessary to sample a 
selection of disordered structures. Furthermore, we may anticipate that the presence of 
charge will have a greater influence in distorting the cluster from its neutral shape than 
it does in the solid-like regime. We are currently investigating the influence of these 
effects. 

We conclude this section by noting that in some systems there may be more 
sophisticated mechanisms that affect the equilibrium distribution of charges. For 
example, the surface of liquid water may develop a slight excess charge even without 
dissolved ions due to the asymmetry of hydrogen bonds, which results in a small transfer 
of charge between molecules [32, 33]. A related effect can lead to some of the charge of 
a solvated ion appearing near the water surface via a chain of partial charge transfers 
through oriented water molecules [7]. Hence, some of the charge may appear at the 
surface even though the ions are located within the bulk liquid. Our simple Lennard- 
Jones model does not attempt to capture the structure of water or the rather specific 
effects associated with hydrogen bonding. However, the model could be modified to 
incorporate some degree of charge separation within individual particles by adding 
a polarisability term to the potential. Self-consistent solution of the polarisability 
equations to obtain the induced dipole moments adds a considerable computational 
overhead to the simulations. Some preliminary calculations of this sort reveal, as 



by 



n >n (s)=q* 2 W N>n (s). 



PN,n( S ) X ^N,n( S ) eX P ~ V N,n( S )/ kT J 



(4) 
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Figure 6. Probability distributions of the number of charged particles lying on the 
cluster surface for various combinations of parameters as listed in the legend. For each 
pair of curves, the solid line is from explicit MC calculations while the dashed line is 
the prediction of Eq. (4). The five sets of data correspond to the following parameters: 



n = 30, q* = 
n = 50, q* = 
n = 100, q* 
n = 100, q* 
n = 150, q* 



0.33, T* = 
0.16, T* = 
= 0.14, T* 
= 0.44, T* 
= 0.336, T 



0.2, Q*/Q* 
0.3, Q I Q ma> 
= 0.2, Q*/Q* mi 
= 0.2, Q*/Q* m: 
= 0.2, Q*/Q* 



0.28 

0.20 
= 0.26 
= 0.80 

= 0.80 



expected, that polarisability leads to an even lower tendency of charged particles to 
lie at the surface of the droplet. 

5. Concluding remarks 

We have taken a thermodynamic approach to the stability and structure of the 309- 
particle Lennard- Jones cluster bearing charge on some of the particles. By confining 
the cluster to a spherical container we have calculated the equilibrium between intact 
and dissociated configurations and have obtained statistics for the occupation of surface 
sites by charged particles. 

The key observation from this work is that both the stability of the cluster and the 
extent to which charges are likely to be found at its surface depend not just on the total 
charge but also on the number of particles over which the charge is distributed. When 
the charge is spread over more particles, the cluster is able to sustain a greater total 
charge but the fraction of the charges that reside at surface sites is lower. In general, 
the cluster is prone to break up at charges lower than those required to drive all the 
charged particles to the surface. 

The location of charged particles in a sphere became an important topic following 
the introduction of the "plum pudding" model of the atom by Thomson in 1904 [34]. 
Since then, the "Thomson problem" has come to refer to finding the energetically 
optimal distribution of point charges constrained to the surface of a sphere. This 
idealised model has received considerable attention as an interesting exercise in global 
optimisation [35]. Although our charged Lennard- Jones cluster seems to share some 
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common features with the Thomson problem, we note that the charged particles would 
only be free to move continuously on the cluster surface when the cluster has melted, 
since at low temperatures the charges would be constrained to lie at the sites of the 
icosahedral structure. However, at liquid-like temperatures we have seen that the 
charges do not all reside at the surface of the sphere and that entropy plays a large 
role in determining their distribution. This means that it is not easy to make contact 
with the Thomson problem in this context. 

Returning to the original motivation of a better understanding of electrospray 
ionisation, it will now be essential to move to a dynamical treatment of cluster instability 
and the mechanisms of fragmentation. An unconfined cluster undergoes evaporation 
even of neutral particles, which would have the effect of concentrating charge within a 
shrinking cluster. It will therefore be important to approach the instability from below 
as well as looking at the highly- charged regime. The Lennard- Jones model should allow 
such processes to be investigated in detail for a cluster bound by simple interactions. 
We intend to develop this coarse-grained approach in a number of directions, including 
inhomogeneous droplets in which a solute has been introduced. 
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